Econometria III

Modelos de Escolha Qualitativa

Rafael Bressan

Modelos de Escolha Qualitativa

Até agora, nossos modelos ficaram assim: \[ \begin{align} y &= b_0 + b_1 x + e \\ e &\sim D\left(0,\sigma^2\right) \end{align} \]

  • A suposição de distribuição em \(e\):

  • Em princípio implica que \(y \in \mathbb{R}\).

  • Resultados de testes, renda, taxas de criminalidade, etc. são todos resultados contínuos. ✅

  • Mas alguns resultados são claramente binários (ou seja, VERDADEIRO ou FALSO):

  • Ou você trabalha ou não,

  • Ou você tem filhos ou não,

  • Ou você comprou um produto ou não,

  • Você jogou uma moeda e saiu cara ou coroa.

Resultados Binários

  • Resultados restritos a FALSO vs VERDADEIRO, ou 0 vs 1.

  • Teríamos \(y \in \{0,1\}\).

  • Nessas situações, estamos principalmente interessados em estimar a probabilidade de resposta ou a probabilidade de sucesso:

\[p(x) = \Pr(y=1 | x)\]

  • como \(p(x)\) muda quando mudamos \(x\)?

    Se aumentarmos \(x\) em uma unidade, como a probabilidade de \(y=1\) mudaria?

Lembrando o Experimento de Bernoulli

Lembre-se da Distribuição de Bernoulli?: Chamamos uma variável aleatória \(y \in \{0,1\}\) tal que

\[ \begin{align} \Pr(y = 1) &= p \\ \Pr(y = 0) &= 1-p \\ p &\in[0,1] \end{align} \]

uma variável aleatória de Bernoulli.

Condicione essas probabilidades em uma covariada \(x\)

\[ \begin{align} \Pr(y = 1 | X = x) &= p(x) \\ \Pr(y = 0 | X = x) &= 1-p(x) \\ p(x) &\in[0,1] \end{align} \]

  • Particularmente: valor esperado (ou seja, a média) de \(Y\) dado \(x\) \[ E[y | x] = 1 \times p(x) + 0 \times (1-p(x)) = p(x) \]

  • Muitas vezes modelamos expectativas condicionais 😉

O Modelo de Probabilidade Linear (MPL)

  • A opção mais simples. Modele a probabilidade de resposta como \[ \Pr(y = 1 | x) = p(x) = \beta_0 + \beta_1 x_1 + \dots + \beta_K x_K \]

  • Interpretação: uma mudança de 1 unidade em \(x_1\), resulta em uma mudança de \(\beta_1\) em \(p(x)\).

Exemplo: Mroz (1987)

  • Participação feminina no mercado de trabalho

  • Como o status de inlf (na força de trabalho) depende da renda familiar da mulher solteira, sua educação, idade e número de filhos pequenos?

Mroz 1987

Code
mroz = woo.dataWoo("mroz")
# Filled histogram. Colors by inlf. Ages in bins of 2 years. Use seaborn
sns.histplot(
    mroz, x="age", hue="inlf", multiple="fill", bins=range(0, 100, 2), stat="density"
)
# Name axes: x = Idade, y = Proporção
plt.xlabel("Idade")
plt.ylabel("Proporção")
Text(0, 0.5, 'Proporção')

Rodando o MPL

Code
lpm = smf.ols(
    "inlf ~ nwifeinc + educ + exper + I(exper**2) + age + I(age**2) + kidslt6",
    data=mroz,
).fit()
df = summary_col([lpm], model_names=["LPM"]).tables[0]
df.iloc[[0, 1, 2, 3, 14, 15, 17], :]
LPM
Intercept 0.3220
(0.4864)
nwifeinc -0.0034
(0.0015)
kidslt6 -0.2604
(0.0341)
R-squared Adj. 0.2568
  • idêntico aos nossos modelos de regressão linear anteriores

  • Variável dependente inlf assume somente dois valores, 0 ou 1.

  • Resultados: se a renda da mulher solteira aumenta em 10 (ou seja, 10.000 USD), \(p(x)\) cai em 0,034 (isso é um efeito pequeno!),

  • uma criança pequena adicional reduziria a probabilidade de trabalho em 0,26 (isso é grande).

  • Até agora, tudo simples.️

MPL: Prevendo probabilidades negativas?!

Code
pr = lpm.predict()
sns.scatterplot(x=np.arange(len(pr)), y=np.sort(pr))
sns.lineplot(x=[0, len(pr)], y=[0, 0], color="red", linewidth=2)
sns.lineplot(x=[0, len(pr)], y=[1, 1], color="red", linewidth=2)
plt.ylabel("P(inlf = 1)")
plt.xlabel("Índice")
plt.title("Previsões ordenadas do MPL")
Text(0.5, 1.0, 'Previsões ordenadas do MPL')



* As previsões do MPL de \(p(x)\) não estão garantidas no intervalo unitário \([0,1]\).

  • Lembre-se: \(e \sim D\left(0,\sigma^2\right)\)

  • aqui, algumas probabilidades menores que zero!

  • Particularmente irritante se você quiser previsões: O que é probabilidade -0,3? 🤔

MPL no modelo saturado: sem problemas!

Code
mroz["age_fct"] = pd.cut(mroz["age"], bins=3, labels=False)
mroz["huswage_fct"] = pd.cut(mroz["huswage"], bins=2, labels=False)
mroz["classes"] = (
    "age_" + mroz["age_fct"].astype(str) + "_hus_" + mroz["huswage_fct"].astype(str)
)

lpm_saturated = smf.ols("inlf ~ -1 + classes", data=mroz).fit()
lpm_saturated.summary(slim=True).tables[1]
coef std err t P>|t| [0.025 0.975]
classes[age_0_hus_0] 0.6108 0.028 22.029 0.000 0.556 0.665
classes[age_0_hus_1] 0 0.349 0 1.000 -0.684 0.684
classes[age_1_hus_0] 0.5851 0.029 19.936 0.000 0.527 0.643
classes[age_1_hus_1] 0.3333 0.201 1.657 0.098 -0.062 0.728
classes[age_2_hus_0] 0.4621 0.041 11.289 0.000 0.382 0.542
classes[age_2_hus_1] 0.5000 0.349 1.435 0.152 -0.184 1.184
  • modelo saturado : só tem variáveis explicativas binárias (dummies)

  • Cada classe: \(p(x)\) dentro daquela célula.

MPL no modelo saturado: sem problemas!

Text(0.5, 0, 'Coeficiente')

  • Cada estimativa pontual: \(p(x)\) dentro do intervalo \([0,1]\).

  • Mulheres da faixa etária mais jovem e de menor renda do marido (coeficiente age_0_hus_0) têm a maior probabilidade de trabalhar (0.611).

Modelos de Resposta Binária Não-Lineares

Nesta classe de modelos mudamos a forma como modelamos a probabilidade de resposta \(p(x)\). Em vez da estrutura linear simples de cima, escrevemos \[ \Pr(y = 1 | x) = p(x) = G \left(\beta_0 + \beta_1 x_1 + \dots + \beta_K x_K \right) \]

  • quase idêntico ao MPL!

  • exceto que o índice linear \(\beta_0 + \beta_1 x_1 + \dots + \beta_K x_K\) agora está dentro da função de ligação \(G(\cdot)\) (i.e. link function).

  • Propriedade principal de \(G\): transforma qualquer \(z\in \mathbb{R}\) em um número no intervalo \((0,1)\).

  • Isso resolve nosso problema de previsões fora do intervalo \(\left[0, 1\right]\) para probabilidades.

Modelos de Resposta Binária Não-Lineares

\(G\): probit e logit


Para probit e logit:

  1. qualquer valor \(x\) resulta em um valor \(p(x)\) entre 0 e 1.

  2. estritamente crescentes.

  3. Logit tem caudas mais longas.

Modelos de Resposta Binária Não-Lineares

\(G\): probit e logit


Probit:

\(G(z)=\Phi(z) = \int_{-\infty}^z \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}t^2}dt\) Distribuição Normal

Logit:

\(G(z)=\Lambda(z)=\frac{1}{1+\exp(-z)}\) Distribuição Logística

Modelos de Resposta Binária Não-Lineares

Rodando probit e logit no Python

  • Podemos usar a função glm para rodar um modelo linear generalizado

  • Isso generaliza nosso modelo linear padrão. Temos que especificar uma família e um link.

  • Porém o statsmodels possui funções especializadas para probit e logit.

Code
probit = smf.probit("inlf ~ age", data=mroz).fit()

logit = smf.logit("inlf ~ age", data=mroz).fit()
Optimization terminated successfully.
         Current function value: 0.680534
         Iterations 4
Optimization terminated successfully.
         Current function value: 0.680519
         Iterations 4

Interpretação

Probit Logit
Intercept 0.7065*** 1.1364***
(0.2470) (0.3983)
age -0.0125** -0.0202**
(0.0057) (0.0092)

Standard errors in parentheses.
* p<.1, ** p<.05, ***p<.01
  • coeficiente de probit para idade é -0.013

  • -0.02 para logit,

  • impacto da idade na probabilidade de trabalhar é negativo

  • No entanto, quão negativo? Não podemos dizer!

Interpretação

\[ \Pr(y = 1 | \text{age})= G \left(x \beta\right) = G \left(\beta_0 + \beta_1 \text{age} \right) \] e o efeito marginal de idade na probabilidade de resposta positiva é

\[\frac{\partial{\Pr(y = 1 | \text{age})}}{ \partial{\text{age}}} = g \left(\beta_0 + \beta_1 \text{age} \right) \beta_1\]

  • \(g\) é definida como \(g(z) = \frac{dG}{dz}(z)\) - a primeira derivada de \(G\) (sendo \(G\) uma distribuição, \(g\) é a função densidade).

  • dado que \(G\) é não-linear, isso significa que \(g\) não será constante. Você pode experimentar isso usando este aplicativo aqui:

Interpretação

Não há um único efeito marginal nesses modelos, pois isso depende de onde avaliamos a expressão anterior. Na prática, existem duas abordagens comuns:

  1. reporte o efeito parcial na média (PEA): \[PEA(X_j)=g(\bar{x} \beta) \beta_j\]

  2. relate o efeito parcial médio (APE): \[APE(X_j)=\frac{1}{n} \sum_{i=1}^N g(x_i \beta) \beta_j\]

Felizmente, a classe dos resultados de um GLM tem um método get_margeff que calcula esses efeitos para nós.

Efeitos Marginais

Code
def get_mfx(model):
    return pd.concat(
        [
            model.get_margeff(at=p).summary_frame().loc[:, ["dy/dx", "Std. Err."]]
            for p in ["mean", "overall"]
        ],
        axis=1,
    )


f = "inlf ~ age + kidslt6 + nwifeinc"  # Regressao estimada
probit_reg = smf.probit(f, data=mroz).fit()
logit_reg = smf.logit(f, data=mroz).fit()

probit_mfx = get_mfx(probit_reg)
logit_mfx = get_mfx(logit_reg)
Optimization terminated successfully.
         Current function value: 0.635318
         Iterations 5
Optimization terminated successfully.
         Current function value: 0.635295
         Iterations 5

Efeitos Marginais

Probit

PEA APE
dy/dx Std. Err. dy/dx Std. Err.
age -0.013671 0.002610 -0.012646 0.002284
kidslt6 -0.313910 0.043171 -0.290372 0.035530
nwifeinc -0.004496 0.001629 -0.004159 0.001484

Logit

PEA APE
dy/dx Std. Err. dy/dx Std. Err.
age -0.013926 0.002686 -0.012652 0.002282
kidslt6 -0.321692 0.046099 -0.292258 0.036455
nwifeinc -0.004593 0.001688 -0.004173 0.001505

Qualidade do Ajuste em Modelos Binários

  • Não existe \(R^2\) universalmente aceito para modelos binários.

  • Podemos pensar em um pseudo \(R^2\) que compara nosso modelo com outro sem regressores:

Code
probit_reg.prsquared
0.07084971547449481
  • Mas isso não é super informativo (ao contrário do \(R^2\)). As alterações no valor da log-verossimilhança são altamente não lineares.

Qualidade do Ajuste em Modelos Binários

  • Vamos verificar precisão - qual é a proporção prevista corretamente!

  • Podemos atribuir a previsão 1 caso \(p(x) > 0.5\) e 0 caso contrário.

Code
confusion_tbl = probit_reg.pred_table(threshold=0.5) / probit_reg.nobs

confusion_matrix = pd.DataFrame(
    confusion_tbl,
    index=["Real 0", "Real 1"],
    columns=["Previsto 0", "Previsto 1"],
)
confusion_matrix
Previsto 0 Previsto 1
Real 0 0.169987 0.261620
Real 1 0.122178 0.446215

Curvas ROC (Receiver Operating Characteristics)

Nomenclatura

Taxa de Verdadeiros Positivos (TPR, sensitividade, recall): \(TP/P=\) 0.785

Taxa de Falsos Positivos (FPR): \(FP/N=\) 0.606

Acurácia: \((TP + TN)/ (P+N)=\) 0.616

  Previsto 0 Previsto 1
Real 0 0.169987 0.261620
Real 1 0.122178 0.446215

Métrica relacionada a Poder

  Previsto 0 Previsto 1
Real 0 0.169987 0.261620
Real 1 0.122178 0.446215

Métrica relacionada a Erro Tipo I

  Previsto 0 Previsto 1
Real 0 0.169987 0.261620
Real 1 0.122178 0.446215

Curvas ROC

  • O corte de 0,5 é arbitrário. E se todas as probabilidades previstas forem \(> 0,5\), mas nos dados houver cerca de 50% de zeros?

  • Vamos escolher um corte arbitrário \(c \in (0,1)\) e verificar a precisão de cada valor. Isso dá uma visão melhor.

  • Podemos confrontar a taxa de verdadeiros positivos (TPR) com a taxa de falsos positivos (FPR), para cada valor de corte \(c\in[0, 1]\).

    1. TPR: número de mulheres corretamente previstas para trabalhar dividido pelo número de mulheres que trabalham.

    2. FPR: número de mulheres incorretamente previstas para trabalhar dividido pelo número de mulheres não trabalhadoras.

  • Plotar FPR vs TPR para cada \(c\) define a curva ROC.

  • Um bom modelo tem uma curva ROC no canto superior esquerdo: FPR = 0, TPR = 1.

Curvas ROC

Code
from sklearn.metrics import roc_curve, auc, accuracy_score
pred_prob = probit_reg.predict()
fpr, tpr, thresholds = roc_curve(mroz["inlf"], pred_prob)
roc_auc = auc(fpr, tpr)
accuracies = [accuracy_score(mroz["inlf"], pred_prob > t) for t in thresholds]
acc_t = thresholds[np.argmax(accuracies)] # threshold que maximiza a acurácia

plt.plot(fpr, tpr, color="darkorange", lw=2, label=f"ROC curve (area = {roc_auc: .2f})")
plt.legend()
plt.ylabel("TPR")
plt.xlabel("FPR")
Text(0.5, 0, 'FPR')



  • Melhor precisão em torno de \(c=0.54\)

  • ROC sempre acima da linha de 45 graus. Melhor do que atribuição aleatória (jogar uma moeda)!

Estimação de Máxima Verossimilhança

  • Ao invés de usar média e variância condicionais, usaremos a distribuição condicional completa

  • Amostra iid \(\lbrace y_i, x_i \rbrace_{i=1}^N\). Queremos estimar a distribuição \(f(y|x)\)

  • Hipótese: esta distribuição é conhecida, a não ser por um número finito de parâmetros fixos.

    • Impomos um modelo paramétrico para a densidade condicional

Estimação de Máxima Verossimilhança

Exemplo Probit

  • \(y_i^*=x_i\theta + \varepsilon_i\), sendo \(\varepsilon_i\sim N(0,1)\)

  • Não observamos \(y_i^*\), apenas \(y_i\)

\[ \begin{equation*} y_i = \begin{cases} 1 \text{ se } y_i^* > 0\\ 0 \text{ se } y_i^* \leq 0 \end{cases} \end{equation*} \]

  • Dadas as especificações acima: \[ \begin{aligned} P(y_i=1|x_i)&=P(y_i* > 0|x_i)\\ &=P(\varepsilon_i > -x_i\theta |x_i)\\ &=1-\Phi(-x_i\theta)\\ &=\Phi(x_i\theta) \end{aligned} \]

Estimação de Máxima Verossimilhança

Exemplo Probit

  • A função de probabilidade condicional pode ser escrita como: \[p(y_i|x_i)=\Phi(x_i\theta)^{y_i}\cdot [1-\Phi(x_i\theta)]^{1-y_i}\]

  • Esta é a probabilidade de ocorrência de um ponto \((x_i, y_i)\). Para o conjunto de dados observados temos a função de verossimilhança \[\ell(\theta)=\Pi_{i=1}^N \Phi(x_i\theta)^{y_i}\cdot [1-\Phi(x_i\theta)]^{1-y_i}\]

  • Máxima Verossimilhança consiste em maximizar a função de verossimilhança em relação aos parâmetros, \(\theta\).

  • Podemos maximizar a log-verossimilhança por praticidade. \(\mathcal{L}(\theta)=\log(\ell(\theta))=\sum_{i=1}^N y_i\log(\Phi(x_i\theta))+(1-y_i)\log(1-\Phi(x_i\theta))\)

Estimação de Máxima Verossimilhança

Exemplo Probit

\[\hat\theta = \arg \max_{\theta} \mathcal{L}(\theta)\]

  • O vetor de score é dado pela derivada da log-verossimilhança em relação a cada um dos parâmetros \[s(\theta)=\nabla_\theta \mathcal{L}(\theta)=\left[\frac{\partial\mathcal{L}(\theta)}{\partial\theta_1}, \ldots, \frac{\partial\mathcal{L}(\theta)}{\partial\theta_k}\right]^{\prime}\]

  • A Hessiana é a matriz de segundas derivadas \[H(\theta)=\nabla^2_\theta \mathcal{L}(\theta)\]

Estimação de Máxima Verossimilhança

Exemplo Probit

Exercício: Calcule o score do modelo probit.

Estimação de Máxima Verossimilhança

Exemplo Probit - Solução

\(\mathcal{L}(\theta)=\sum_{i=1}^N y_i\log(\Phi(x_i\theta))+(1-y_i)\log(1-\Phi(x_i\theta))\)

\[\begin{align} s(\theta)&=\frac{d\mathcal{L}(\theta)}{d\theta}=\sum_{i=1}^N y_i x_i\frac{\phi(x_i \theta)}{\Phi(x_i \theta)} - (1-y_i)x_i\frac{\phi(x_i \theta)}{1-\Phi(x_i \theta)}\\ &=\sum_{i=1}^N \frac{[1-\Phi(x_i \theta)]y_i x_i\phi(x_i \theta) - \Phi(x_i \theta)(1-y_i)x_i\phi(x_i \theta)}{\Phi(x_i \theta) [1-\Phi(x_i \theta)]}\\ &=\sum_{i=1}^N \frac{x_i\phi(x_i \theta)}{\Phi(x_i \theta) [1-\Phi(x_i \theta)]}[y_i-\Phi(x_i \theta)] \end{align}\]

Estimação de Máxima Verossimilhança

Exemplo Exponencial

  • Suponha que tenhamos um conjunto simples de dados \(\{x_i\}_{i=1}^N\) \(iid\) e queremos ajustar uma distribuição exponencial

  • \(x_i\sim Exp(\lambda)\). \(f(x; \lambda)=\frac{1}{\lambda}e^{-x_i/\lambda}\). \(x_i\geq 0\).

  • Qual é a função de log-verossimilhança?

  • Encontre o estimador de máxima verossimilhança \(\hat\lambda_{mle}\)

Estimação de Máxima Verossimilhança

Exemplo Exponencial

  • Função de verossimilhança: \(L(\lambda)=\Pi_{i=1}^N f(x_i; \lambda)=\Pi_{i=1}^N \frac{1}{\lambda}e^{-x_i/\lambda}\)

  • Log-verossimilhança: \(\mathcal{L}(\lambda)=\sum_{i=1}^N \log \frac{1}{\lambda}e^{-x_i/\lambda}\)

  • Ao final encontramos …

\[\hat\lambda_{mle}=\frac{1}{N}\sum_{i=1}^N x_i\]

Método dos Momentos Generalizado

Método dos Momentos (MM)

  • Podemos utilizar momentos populacionais para identificar parâmetros

    • \(E[u_i]=0\) é um exemplo
  • De fato, tanto MQO quanto VI podem ser entendidos como estimadores MM!

  • Quais são os momentos utilizados no MQO?

  • \(E[u]=0\) e \(E[xu]=0\), que são as equações normais do MQO

  • O MM resolve o análogo amostral destas equações

Método dos Momentos Generalizado

Método dos Momentos (MM)

  • Vamos considerar um setup mais geral:
  • \(\theta_0\in\mathbb{R}^p\) é um vetor dos parâmetros da população
  • \(m(\theta_0)\equiv E[g(Z_i, \theta_0)]=0\), são as equações de momento
  • \(g(Z_i, \theta)\) é uma função vetorial \(r\times 1\)
  • \(Z_i\equiv (X_i, Y_i)\) é uma observação do conjunto de dados
  • \(\theta_0\) é a única solução para a equação de momento
  • MM vai resolver o análogo amostral:

\[\hat{m}(\hat\theta)=\frac{1}{N}\sum_{i=1}^N g(Z_i, \hat\theta)=0\]

Método dos Momentos Generalizado

Clarificando

  • Considere o MQO. Podemos escrever duas equação de momento que identificam os parâmetros do modelo \[ \hat{m}(\hat\beta)=\begin{cases} \frac{1}{N}\sum_{i=1}^N(y_i-\hat\beta_0-\hat\beta_1 x_i)=0\\ \frac{1}{N}\sum_{i=1}^Nx_i(y_i-\hat\beta_0-\hat\beta_1 x_i)=0 \end{cases} \]

  • O MM escreve momentos populacionais em termos dos parâmetros e resolve o sistema de equações, simples!

Método dos Momentos Generalizado

GMM

  • Podemos generalizar o método dos momentos. Daí o nome generalized method of moments - GMM

  • Minimizar a norma \({\displaystyle \|{\hat {m}}(\theta )\|_{W}^{2}={\hat {m}}(\theta )^{\mathsf {T}}\,W{\hat {m}}(\theta)}\)

  • Onde \(W\) é uma matriz positiva definida

  • \({\displaystyle {\hat {\theta }}=\operatorname {arg} \min _{\theta \in \Theta }{\bigg (}{\frac {1}{N}}\sum _{i=1}^{N}g(Z_{i},\theta ){\bigg )}^{\mathsf {T}}{W}{\bigg (}{\frac {1}{N}}\sum_{i=1}^{N}g(Z_{i},\theta ){\bigg )}}\)

  • Especialmente útil quando temos mais equações de momento do que parâmetros

Método dos Momentos Generalizado

Aplicações de GMM

  • MQO, VI e MLE podem ser obtidos a partir de um GMM

  • Aplicações diversas em

    • Cross-Section: modelos não-lineares com variáveis endógenas
    • Séries Temporais: correlação serial e heterocedasticidade
    • Painel: extensões do modelo de efeitos fixos
  • Veja mais em Wooldridge (2001) JEP

📚 Leitura Recomendada

  • WOOLDRIDGE, Jeffrey M. Introdução à econometria: uma abordagem moderna. São Paulo: Cengage Learning, 2016. Tradução da 4ª edição norte-americana por José Antonio Ferreira. Capítulo 17 Modelo com Variáveis Dependentes Limitadas e Correções da Seleção Amostral.

  • GUJARATI, Damodar N.; PORTER, Dawn C. Econometria básica. Porto Alegre: Amgh Editora, 2011. - 5. ed. Capítulo 15 Modelos de regressão de resposta qualitativa

  • WOOLDRIDGE, Jeffrey M. Econometric Analysis of Cross Section and Panel Data. MIT press, 2010. Second Edition. Chapter 15 Binary Response Models

ATÉ A PRÓXIMA AULA!